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ABSTRACT 

We develop a technique to measure radial velocities of stars from spectra that present four sets 
of lines. The algorithm is an extension of the two-dimensional cross-correlation method TODCOR 
to four dimensions. It computes the correlation of the observed spectrum against a combination of 
four templates with all possible shifts, and allows also for the derivation of the light ratios of the 
components. After testing the algorithm and demonstrating its ability to measure Doppler shifts 
accurately even under conditions of heavy line blending, we apply it to the case of the quadruple- 
lined system HD 110555. The primary and secondary components of this previously known visual 
binary (p ~ 0'.'4) are each shown to be double-lined spectroscopic binaries with periods of 57 days 
and 76 days, respectively, making the system a hierarchical quadruple. The secondary in the 76-day 
subsystem contributes only 2.5% to the total light, illustrating the ability of the method to measure 
velocities of very faint components. 

Subject headings: binaries: general — binaries: spectroscopic — binaries: visual — methods: data 
analysis — stars: individual (HD 110555) — techniques: spectroscopic 



1. INTRODUCTION 

The use of digital cross-correlation as an analysis 
tool in spectroscopy dat es back at leas t three d ecades . 
Early applicat i ons b y ISimkinl (| 19741) . lLacvl (|1977T ). 
iTonrv fe Davisl (|1979f ). and others, showed the power 
of the method for measuring accurate radial velocities, 
and with many improvements and refinements over the 
years the technique enjoys widespread use today (for fur- 
ther details and a histori cal perspective see, e.g., iHilll 
119931 iKurtz fe Minkl 119981 ). Compared to the classical 
method of measuring the positions of individual spec- 
tral lines (still occasionally used) , cross-correlation offers 
a number of important advantages beyond expediency, 
principally its ability to make efficient use of all the in- 
formation available in the spectral window. It is thus 
ideal for the analysis of low signal-to- noise spectra, where 
classical methods tend to give poorer results. It is also 
commonly used on composite spectra, e.g., double-lined 
spectroscopic binaries, where the Doppler shifts are typi- 
cally determined from the centroids of the two correlation 
peaks that correspond to the each of the components of 
the binary. Difficulties arise, however, when the peaks 
are not well separated, and this can lead to systematic 
errors in the radial velocities that bias the amplitudes of 
the velocity curves and any physical quantities derived 
from them (such as the masses), or it can prevent the 
measurement altogether. 

To address this problem I Zucker fc Mazehl (| 19941 ) devel- 
oped an extension of the standard one-dimensional cross- 
correlation technique to two dimensions (TODCOR), 
taking advantage of the properties of the Fourier trans- 
form. TODCOR uses two templates, one for each com- 
ponent of the binary, and thus the cross-correlation func- 
tion (CCF) now depends on the velocities of both stars 
relative to their templates. The location of the maxi- 
mum in two-dimensional velocity space corresponds to 
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the Doppler shifts of the two components. This effec- 
tively decouples the primary from the secondary so that, 
if the templates are a good match to the stars, blending 
between them is no longer a concern. A further exten- 
sion of TODC OR to three d i mensi ons was subsequently 
introduced by IZucker et al.l (|1995f ) to deal with triple- 
lined spectra, in which line blending is usually more com- 
mon. In this case the standard one-dimensional CCF 
would in principle reveal three peaks (at favorable or- 
bital phases), but the difficulties of measuring their cen- 
troids by classical means can be even greater than in the 
double-lined case. The success of the three-dimensional 
version of TODCOR to overcome these difficulties is il- 
lustrated by a number of applications over the last decade 
since its introduct i on (ITorres et al.lll995tlJha et al. 2000 ; 
Mazeh et ai]|200lUCovino et aUl200lHTorres et alJhooa . 
20061) . 



A recent study by iTokovinin et al.l (|2006f ) has shown 
that spectroscopic binaries of solar-type are very often 
accompanied by more distant components. The fre- 
quency of such multiple systems appears to be very high: 
up to 97% of spectroscopic binaries with periods shorter 
than 3 days show at least one additional companion, ta- 
pering off to 34% for periods lon ger than 12 days (see 
Fig. 14 bv ITokovinin et alj|2006l ). Out of the 161 sys- 
tems in their sample, 64 were found to be triple and 11 
were found to be quadruple. Thus the frequency of sys- 
tems with four stars is not negligible, and it is to be 
expected that some fraction of them will show lines of all 
four stars in the spectra, although they may not always 
be easy to recognize. These will typically be hierarchi- 
cal quadruple systems composed of two relatively tight 
binaries in a wide orbit around each other, with an angu- 
lar separation of order an arc second or less so that they 
are unresolved at the spectrograph slit (or optical fiber 
entrance). The ability to measure the radial velocities 
of the 4 stars in these systems is essential for deriving 
properties such as their masses, as well as other orbital 
characteristics, and these in turn are important for un- 
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derstanding the origin and evolution of multiple systems 
in general, and higher multiplicity systems in particular. 

The extensive spectroscopic surveys at the Harvard- 
Smithsonian Center for Astrophysics ( CfA) with the CfA 
Digital Speedometers (|Lathaml Il992f ) have revealed an 
appreciable number of triple-lined systems, and indeed 
even some cases with four sets of lines that have until now 
remained unsolved because of the complexity of the anal- 
ysis and our concerns about blending. This provides the 
motivation for the present work, in which we develop an 
extension of TODCOR to four dimensions, allowing the 
treatment of such cases. We demonstrate the effective- 
ness of the method with numerical simulations, showing 
that it is possible to measure the velocities of the four 
stars even in situations with severe line blending that 
would be hopeless using standard one-dimensional cross- 
correlation techniques. We then apply the method to 
the real case of HD 110555, which enabled us to detect 
for the first time the very faint 4th component and to 
measure its Doppler shift, in addition to obtaining the 
relative brightnesses of the stars. 

2. DESCRIPTION OF THE METHOD 

TODCOR as introduced by iZucker fc Mazehl (|1994D 
is a two-dimensional cross-correlation scheme in which 
the observed spectra, f(n), are cross-correlated against 
a composite template, g(n), that is the sum of two sep- 
arate templates selected to match the properties of each 
of the binary components. The template g for this two- 
dimensional case is given by g = g\(n — si)+a g2{n — S2), 
where s\ and S2 represent the shifts of the individual 
templates (<?i and (72) needed to match the true Doppler 
shift of the respective component. The coefficient a is the 
light ratio between the secondary and the primary, which 
is assumed initially to be known. The cross-correlation 
between / and g is performed for all possible shifts of 
the two templates, and the resulting function of s\ and 
s 2 will typically have a global maximum, the location of 
which provides the Doppler shifts of both stars. Brute- 
force computation of these cross-correla tions over all pos- 
sible shifts is impractical. However, IZucker fc Mazehl 
( 1994) showed that the two-dimensional correlation func- 
tion , S2) can be expressed analytically 
in terms of three one-dimensional correlation functions 
which are much more economical to compute. As they 
pointed out, computing one-dimensional CCFs is an 
0(N 2 ) process, where N is the number of pixels in /, 
but by using the properties of Fourier transforms the 
problem can be reduced to an O(NlogN) process. This 
reduction in the number of operations required is pre- 
served in TZ^ 2 \ and is what makes the method practi- 
cal. In the more general case in which the light ratio 
is unknown, the correlation function 1ZS 2 ' also depends 
on a, or 1Z^ — TZ^ 2 \sx,S2 1 ct)- An analytical expres- 
sion can be found for a that maximizes the correlation 
for each pair of shifts s\ and S2, and this expression 
depends on the same three one-dimensional CCFs (see 
IZucker fc Mazehl Il99l . 

The ext ension of TODCOR to three dimensions pre- 
sented by IZucker et alj (|1995f ) is fairly straightforward, 
the algebra being somewhat more involved. The Fourier 
transform properties once again reduce the computa- 
tional problem to an 0(N log N) process. Conceptually 
it is just as straightforward to extend the scheme to four 



dimensions for the analysis of quadruple-lined spectro- 
scopic systems, although the algebra in this case is con- 
siderably more complex. The template is now assumed 
to be a combination of four separate (i.e., possibly differ- 
ent) templates, each with its own shift, and three light 
ratios: 

g = 9i{n-si) + a g2(n~s 2 )+P ff3(n-s 3 )+7 g A (n~s A ) , 

where (3 and 7 are the light ratio of the tertiary and qua- 
ternary components relative to the primary. The four- 
dimensional CCF for the case in which the light ratios 
are known, 7Z^ = 1Z^ (si, S2, S3, S4), can again be ex- 
pressed analytically in terms of one-dimensional correla- 
tion functions and three light ratios, as described in the 
Appendix. These one-dimensional correlations are be- 
tween the observed spectrum and each of the four tem- 
plates, as well as between pairs of templates (in all com- 
binations). To determine the velocities of the four stars 
one simply searches for the maximum of IZ^ in the 
four-dimensional space of {si, S2, S3, S4}. For the more 
likely case in which the light ratios are unknown to be- 
gin with, analytical expressions can be found for a, (3, 
and 7 in terms of the same one-dimensional CCFs men- 
tioned above (see Appendix). These values can then be 
inserted into the expression for"K (4) . Thus the method 
is completely general in that it permits in principle the 
determination of all unknown quantities directly from 
the observed spectra. The only requirement is that the 
four templates chosen for the cross-correlations be a rea- 
sonably good match to each of the components of the 
quadruple system. 

3. TESTING THE ALGORITHM 

In order to test the method we generated syn- 
thetic composite spectra to simulate observations of a 
quadruple-lined system by adding together four calcu- 
lated spectra with various relative shifts and intensity 
ratios. Those same four calculated spectra were in turn 
used as individual templates (.91, <?2, .93, 34) in our four- 
dimensional extension of TODCOR. All synthetic spec- 
tra were taken from a large libra ry of calculated spectra 
created by Jon Mo rse (see also iNordstrom et al.l 11994 
lLatham et al.ll2002T ). based on model atmospheres by R. 
L. Kurucz 1 , which we also use in the next section to ana- 
lyze a real case. These spectra cover a wavelength range 
of approximately 80 A centered at 5187 A, and have a 
resolving power of A/AA « 35,000. They were intended 
for use with spectroscopic observations obtained with the 
CfA Digital Speedometers, which have a narrower wave- 
length coverage of only 45 A, so for comparison purposes 
we have considered only this reduced spectral window in 
our tests. Furthermore, in order to make the simulations 
more realistic we have added Poisson noise corresponding 
to different signal-to- noise ratios (SNR). 

As a first example among the many experiments we 
carried out, we consider a case with the relative shifts of 
the four individual stars selected to demonstrate the abil- 
ity of the algorithm to measure velocities under condi- 
tions of significant line blending. Our simulated quadru- 
ple system in this case is composed of non-rotating stars 
of spectral types GOV, G2V, GOV, and G2V for the pri- 
mary, secondary, tertiary, and quaternary, respectively, 

1 Available at http://cfaku5.cfa.harvard.edu 
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corresponding to effective temperatures (T c ff) for the 
templates of 6000 K, 5750 K, 6000 K, and 5750 K. 
The templates were added together (initially with no 
noise) with Doppler shifts of +20 kms -1 , —20 kms -1 , 
+ 10 kms -1 , and —5 kms -1 , respectively, and with in- 
tensity ratios of a = 0.50 (secondary/primary), j3 = 0.80 
(tertiary/primary), and 7 = 0.50 (quaternary/primary). 
Although some signs of the composite nature of this arti- 
ficial spectrum can be seen by careful visual comparison 
with the primary template, it is not possible to identify 
the lines of the four components, and this would only 
worsen in the presence of noise. The top panel of Figure[T] 
shows the standard one-dimensional CCF we obtain by 
correlating the above noiseless synthetic composite spec- 
trum against the primary template. This and all other 
one-dimensional CCFs in t his paper are comput ed us- 
ing the IRAF 2 task XCSA0 (jKurtz k. Minkl I1998D . The 
input velocities of each component are indicated in the 
figure with vertical dotted lines. The blending of the cor- 
relation peaks is such that the velocities of the individ- 
ual stars are not measurable in this case with standard 
techniques. Application of the four-dimensional cross- 
correlation technique, on the other hand, is able to re- 
cover the four velocities as well as the light ratios reli- 
ably. To show this quantitatively and to provide at the 
same time a realistic measure of the uncertainty due to 
noise, we have added Poisson noise to the input compos- 
ite spectrum corresponding to a SNR = 25, and gener- 
ated 50 quadruple-lined spectra with different noise but 
with the same shifts. We then processed these spectra 
as described above. The results are shown in Table Q] 
(middle section under "Test 1" ) , where we list the mean 
velocities obtained for each star from the 50 spectra, as 
well as the mean light ratios. The uncertainties attached 
to these means are simply the scatter of the 50 values, so 
they represent the typical error for a determination from 
a single spectrum. Both the input velocities and the in- 
put light ratios are recovered to well within their errors, 
showing the power of the technique in this very demand- 
ing case. Cross-sections of the four-dimensional CCF for 
the case of no noise are shown in the lower panels of Fig- 
ure [T] For example, the second panel shows a cut along 
the velocity axis corresponding to the primary star, with 
the velocities of the other three components held fixed 
at the values that maximize the correlation function. A 
peak is clearly seen precisely at the location of the input 
value for the velocity of that star (dotted line). Similarly, 
cross-sections for each of the other components show a 
prominent peak at the correct velocity. 

As a second test we reduced the SNR to 10 and re- 
peated the velocity and light ratio determinations for 
the 50 simulated spectra (see Table [TJ "Test 2"). In 
this case the velocities are typically still recovered quite 
reliably (within a few kms -1 ), and although the mean 
values of a, /3, and 7 are also not far from their input 
values, their scatter is much larger. However, as stated 
above this scatter represents the typical error for a single 
measurement. In a real application there will in general 
be multiple spectra of the object, so if the light ratios 

2 IRAF is distributed by the National Optical Astronomy Ob- 
servatories, which is operated by the Association of Universities 
for Research in Astronomy, Inc., under contract with the National 
Science Foundation. 



are constant one can hope for a y/~N gain in taking the 
averages. 

A third experiment was carried out in which we exam- 
ined the ability of the algorithm to recover faint com- 
ponents in quadruple-lined spectra. For this test we 
adopted templates corresponding to spectral types GOV, 
G2V, G2V, and K4V (temperatures of 6000 K, 5750 K, 
5750 K, and 4750 K), and we relaxed the line blending by 
adopting somewhat larger Doppler shifts of +50 kms" 1 , 
-50 kms' 1 , +20 kms" 1 , and -20 kms" 1 . The SNR 
was set to 50 in this case, and we chose a — 0.70 and 
= 0.60. For 7 we tested a range of values from 0.20 
to 0.01, and we found that the correct velocities were al- 
ways recovered reliably by the algorithm for the primary, 
secondary, and tertiary components, and that the input 
values of a and (3 were recovered as well. The veloci- 
ties for the fourth star deteriorated with decreasing 7, 
as expected, but were still measurable down to 7 = 0.03 
(s 4 = -22.4+1.9 kms" 1 for the mean of 50 trials). At 
this light ratio the recovered value of 7 was 0.039 + 0.004, 
which is slightly overestimated. At a lower light ratio of 
2% the fourth star was only detectable in about half of 
these artificial spectra. 

The above tests and many others we carried out show 
that the algorithm performs well under conditions of se- 
vere line blending allowing the reliable determination of 
radial velocities, and also demonstrate its ability to de- 
tect faint components in quadruple-lined spectra. It is 
clear that the sensitivity to faint companions depends 
strongly on the degree of line blending of the component 
in question, as well as on the SNR of the spectra. In the 
next section we apply the method to a real case. The 
performance of the algorithm will usually depend also on 
how closely the four templates represent the individual 
stars. In the experiments in this section the templates 
we have used are, by construction, a perfect match to 
the components. The effects of mismatch have not been 
addressed here because they depend on the properties 
of the particular system under investigation, and it is 
difficult to make general statements that are useful in a 
different case. We discuss this further in SJH 

4. APPLICATION TO HD 110555 

The relatively bright star HD 110555 (also known 
as BD+06 2647, HIP 62034, and ADS 8639 AB, a = 
12 h 42 m 55!36, 8 = +05°15'59'.'l, J2000, spectral type 
G5, V = 8.36) is a close (p ~ 0'.'4) visual binary dis- 
covere d by R. G. Aitken at the Lick Observatory in 
1907 (IAitkenlfT908h . It was observed spectroscopically 
at the CfA as part of a large sample of G dwarfs that 
have been monitored for more than 15 years. In the 
course of that project we obtained a total of 48 usable 
spectra of HD 110555 between 1993 January and 2004 
March, mostly with an echelle spectrograph mounted on 
the 1.5m Wyeth reflector at the Oak Ridge Observatory 
(Harvard, Massachusetts). A single echelle order span- 
ning 45 A was recorded with a photon-counting Reticon 
diode array at a mean wavelength of 5188.5 A, which 
includes the Mg I b triplet. The SNRs of these observa- 
tions range from about 11 to 42 per resolution element 
of 8.5 kms -1 . Two of the spectra were obtained with a 
nearly identical setup using the 1.5m Tillinghast reflec- 
tor at the F. L. Whipple Observatory (Mount Hopkins, 
Arizona) . Because of the close separation of the binary 
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and the 1" width of the spectrograph slit, the system 
was observed as a single object. The zero-point of the 
velocity scale was monitored by means of exposures of 
the dusk and dawn sky, and small run-to-run corrections 
were applied to t he velocities derived below in the man- 
ner described bv lStefanik et al.l (|1999j ). 

The first few spectra of HD 110555 showed obvious 
signs of at least two sets of lines separated in velocity 
by as much as 80 kms -1 , which, in view of the small 
angular separation of the visual pair, suggested that at 
least one of the visual components is a spectroscopic bi- 
nary with a relatively short period. Figure [2] includes 
a small sampling of our observations showing the spec- 
tra and one-dimensional CCFs computed using a tem- 
plate corresponding to a GO star. Some of the obser- 
vations even showed three sets of lines, and the appear- 
ance of the correlation functions varied on a timescale 
of a week or so. Changes in velocity had not previ- 
ously been measured for HD 110555, and in fact we 
are only aware of a si ngle velocity mea s ureme nt in the 
literature reported by [Nordstro m et al.l (|2004h yielding 
RV = —40.9 ±0.8 kms -1 , which is, however, unlikely to 
be very meaningful in view of the complicated nature of 
the system. 

Preliminary velocity determinations were carried out 
initially using the two-dimensional version of TODCOR, 
focusing on the two dominant correlation peaks, and later 
with the three-dimensional version of the algorithm to ac- 
count for the third peak. For these analyses we adopted 
solar-type templates for all stars (Toff = 5750 K), with 
no rotational broadening since the composite spectra do 
not exhibit significant line broadening at phases where 
the lines are all blended 3 . The surface gravity was set to 
log<? = 4.5, as appropriate for dwarfs, and solar compo- 
sition was assumed. The results showed that the system 
is composed of a double-lined binary with a period of 
about 57 days, and a single-lined binary with a period of 
~76 days. No sign of the secondary of the latter system 
was apparent. The preliminary center-of-mass velocities 
of these orbits were within a few km s _1 of each other, in- 
dicating the physical association of the two binaries and 
leading us to conclude that they most likely correspond 
to each of the components of the visual binary. The 
HD 110555 system is therefore a hierarchical quadruple. 

We then applied our new algorithm in an attempt to 
detect the fourth star (secondary of the 76-day binary) 
and measure its velocity. The template adopted in this 
case was somewhat cooler (T e g = 5000 K) . All 48 of our 
spectra showed evidence of a weak set of lines at approx- 
imately the expected location. An example is shown in 
Figure [H which corresponds to the same spectrum dis- 
played in the bottom panel of Figure [H in which the three 
brighter components are quite heavily blended. The top 
panel of Figure [3] shows the one-dimensional CCF once 
again, only on an expanded scale, and the lower four pan- 
els show cuts of the four-dimensional CCF along the axes 
corresponding to each of the components, as in Figure [TJ 
We indicate with vertical dotted lines the predicted ve- 
locities from our final spectroscopic orbits, described be- 

3 A v sin i measurement of 9 kms -1 was reported for HD 110555 
bv lNordstrom et~a l. (2004), but as with the velocity estimate from 
this source quoted earlier, its interpretation is difficult and we con- 
sider it at best only as an upper limit. 



low. The cross-section for the quaternary is noisy, but 
it does present a peak at the expected location, and the 
same is true for each of our other spectra. 

The preliminary light ratios we derived clearly indi- 
cated that the 57-day binary is brighter than the 76-day 
binary, so it must correspond to the visual primary. We 
refer to this system as "A" , composed of stars Aa and 
Ab, following the usual spectroscopic notation. Similar 
designations are adopted for the visual secondary, "B". 
In order to fine-tune the templates for the four stars we 
made use of model i sochrones from the Padova series by 
iGirardi et alJ ([20001 ) . combined with all available obser- 
vational constraints. These include, in addition to the 
three preliminary light ratios, the combined B — V color 
of the q uadruple sy stem as listed in the Hipparcos Cat- 
alogue ([ESAl Il997^ . the V-J, V-H, and V-K in- 
dices derived from the Jo hnson V magnitude ( Hippar- 
cos) and 2MASS Catalog (jSkrutskie et all 12009) . prop- 
erly converted to the same pho tometric system as the 
stellar models ([Carpenter! 12001). the magnitude differ- 
ence of the visual pair (AV) as measured by Hipparcos 
and c onverted to the visual band following iHarmariecl 
(1998), and the preliminary mass ratios for the two spec- 
troscopic binaries. By requiring simultaneous agreement 
with all observational quantities we were able to select 
four stars from a representative model isochrone of solar 
metallicity and age 3 Gyr that yield a consistent pic- 
ture of the system, and allow us to read off the effective 
temperatures as well as other theoretical properties. The 
temperatures were close to 6000 K for Aa and Ab, 5500 K 
for Ba, and 4750 K for the fainter component Bb, which 
are the nearest values in our library of synthetic spec- 
tra. The surface gravities were not far from logg = 4.5, 
as expected for dwarfs, supporting our earlier choice of 
this value. The fairly close agreement between all ob- 
servational constraints and the results from the model is 
shown in Table [2j where the differences in the last col- 
umn are seen to be of order 1.3er or less 4 . The inferred 
properties of the four stars are listed in Table [3] They 
are unevolved dwarfs, so the stellar characteristics are 
quite insensitive to age for our purposes. Their location 
in the H-R diagram is illustrated in Figure [4] Templates 
with the above parameters and zero rotational broaden- 
ing were then used to derive improved radial velocities, 
and the light ratios and mass ratios also changed slightly. 
One more iteration was sufficient to reach convergence on 
the temperature determination for the templates, given 
the relatively coarse 250 K step in our library of syn- 
thetic spectra. The final velocities were derived with the 
above template parameters, and the light ratios obtained 
were a = 0.92, (3 = 0.38, and 7 = 0.06, determined as 
described in the Appendix. Uncertainties for these ratios 
are estimated to be about 0.03. 

We present the velocities for the visual primary and vi- 
sual secondary components in Table |4] and Table [U The 
uncertainties (a) associated with these measurements re- 
flect both the SNR of each spectrum and the relative 
brightness of each component. They were derived as a 
byproduct of the orbital solutions, by assigning initial 

4 The somewhat larger value of 7 predicted by the isochrone may 
in fact be due at least in part to deficiencies in the models for lower 
mass stars, such as missing opacities in th e optical bands (see, e.g., 
IDelfosse et aL|[2000l; IChabrier et~aTll2005l) . 
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weights to the observations proportional to the strength 
of each exposure, and then scaling these initial errors 
so that the reduced x 2 values are near unity separately 
for the primary and secondary in each subsystem. The 
orbital elements we obtain for the two spectroscopic bi- 
naries are listed in Table [6] They are the period (P), 
center-of-mass velocity (Vb), velocity semi-amplitudes 
(-Kprimj K sec ), eccentricity (e), longitude of periastron 
for the primary (oj pr i m ), and the time of periastron pas- 
sage (T). Derived quantities of interest are listed as 
well (minimum masses, mass ratio, and projected semi- 
major axes). The velocity measurements in these fits 
have been weighted in the usual way according to their 
uncertainties. The rms residuals of the fit for star Ba 
(~1.7 kms -1 ) and especially Bb (~5.4 kms^ 1 ) are con- 
siderably larger than those of Aa and Ab (~0.7 kms -1 ) 
on account of their faintness, which represents only about 
16% and 2.5% of the total brightness of the system at this 
wavelength. Stars Aa and Ab contribute 42.5% and 39% 
to the total light. The velocity measurements and corre- 
sponding orbital fits are shown graphically in Figure 

A visual orbit determi nation for H D 110555 A and 
B has been published by iLingj (|2004f ) with a period of 
1263±125 yr, an angular semimajor axis of l'.'915±0'/022 
(corresponding to ~150 AU), and a very large eccentric- 
ity of e = 0.983 ± 0.001. However, due to the small 
number of measurements 5 and their coverage of only a 
fraction of the orbit, it is assigned "Grade 5" ("indeter- 
minate") in the Sixth Catalog of Orbits of Visual Binar y 
Stars maintained at the USNO (jHartkopf et all 1200 If ), 
and must b e regarded as very preliminary. Based on 
this solution [Ling ( 2004) reported a dynamical mass for 
the system of 2.28 M©, which is, however, smaller than 
the sum of the minimum masses of the four stars from 
our spectroscopic orbits (3.44 ± 0.07 M ; Table®). The 
total mass inferred from our modeling is 3.8 M Q (Ta- 
ble [3|) . It is quite possible that the visual elements can 
in fact be improved by imposing our total mass as a con- 
straint. A further constraint is provided by our radial 
velocities of A and B. They span only 11 yr, and we 
see no trend in the residuals that would reflect motion 
in the outer orbit, but the difference in the center-of- 
mass velocities of the two spectroscopic binary orbits, 
V (B) - V {A) = -0.26 ± 0.23 kms" 1 (Table©, pro- 
vides information that is orthogonal to the astromctry, 
and is thus potentially very important. We note, finally, 
that as a result of our isochrone fitting described above, 
the inclination angles of the two spectroscopic binaries 
are inferred to be ~72° (visual primary) and ^81° (vi- 
sual secondary), with uncertainties of roughly 5°. The 
inclination reported for the visual orbit is rather similar 
(85°). 

5. DISCUSSION 

The cross-correlation technique described and tested in 
this paper, which is an extension of the two-dimensional 
algorithm TODCOR, provides a conceptually simple and 
elegant way of deriving reliable radial velocities from stel- 
lar spectra in which four sets of lines are present. It 
may be thought of as a prescription for cross-correlating 

5 Nine micrometer measurements, one speckle observation, and 
the measurement from the Hipparcos mission. The time span of 
these data is 92 yr. 



the observed spectra against a composite template made 
of the combination of four individual templates (one for 
each component) with all possible shifts. The ability 
to make use of four templates selected to match the 
properties of each star is one of the great strengths of 
the method. Standard one-dimensional cross-correlation 
methods use a single template, which may be a good 
representation of one of the stars but will in general not 
match the others as well. Thus the correlation will never 
be optimal. Another strength of our four-dimensional 
algorithm is its ability to produce good results even 
with low SNR spectra, as illustrated by the example of 
HD 110555, in which the SNRs are as low as 11. Even in 
weak spectra such as these we were able to measure the 
velocity of the faint fourth component that contributes 
only 2.5% to the light of the system (6% of the light of the 
primary), despite seeing no sign of it either in the original 
spectra or in the standard one-dimensional CCFs. The 
method has shown to be very effective as well in dealing 
with line blending, yielding accurate radial velocities for 
HD 110555 although the line separation between two or 
more of the components is sometimes as small as a few 
kms . 

Only a handful of examples of quadruple-lined spec- 
troscopic systems have appeared previously in the lit- 
erature in which reliable radial velocities were de- 
rived for all four star s. To our k nowledge these are 
XY L eo (|Bardenlll987l ). HD 221264 (IWillmitch fe Fekell 
1990j) . ET Boo, VW LMi, and TV UM i lPribulla et ail 
2006li. AO V e l (IGonzalez et all l2006h. and HD 30869 
( Tom kin et alj l2007f ). all of which contain an eclipsing 
pair except for the second and last cases. A number 
of ot her spectroscopic quadrup les have been identified 
(e.g.. lPribulla fc Rucinsk i 2006), although velocity mea- 
surements have not yet been published and it is unclear 
whether these systems show four sets of lines. The ma- 
terial available to the authors of the studies mentioned 
above was generally of much higher quality than that 
employed here, with SNRs typically in the range of 100- 
150. The techniques employed to derive the velocities 
have varied considerably. In the case of XY Leo and 
HD 30869 the centroids of absorption lines were mea- 
sured directly, or the observed spectra were fitted in a 
X 2 sense by coadding spectra of comparison stars ap- 
propriately shifted in velocity, rotationally broadened, 
and scaled in flux. In HD 221264 the lines of the four 
components were separated enough that they could be 
isolated, and radial velocities were derived by focusing 
on a set of unblended lines of one star at a time us- 
ing one-dimensional cross-correlation methods restricted 
to the appropriate wavelength intervals. For ET Boo, 
VW LMi, and TV UMi the authors used t he broad- 
ening function (BF) formalism developed by iRucinskl 
(2002), which reduces the severity of blending compared 
to the one-dimensional cross-correlation technique. The 
peaks in the BFs for these systems were fitted using ei- 
ther Gaussian profiles or rotational profiles. Finally, for 
AO Vel an iterative scheme was used in which the dis- 
entangled spectra of the component s and their Doppler 
shifts were computed in steps (see IGonzalez fc Levatol 
2006). All of these obviously represent useful alterna- 
tives to the technique developed in this paper. 

The advantages or disadvantages of each these methods 
generally depend on the particular case under consider- 
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ation, and their ability to derive accurate radial veloci- 
ties is a strong function of the degree of line blending, 
the relative brightness of the components, and the qual- 
ity of the spectroscopic material, among other factors. 
As demonstrated with the experiments in $3] the four- 
dimensional extension of TODCOR presented here fares 
very well even when the spectra have relatively low SNR. 
Other procedures that have been developed more specif- 
ically for disentangling the spectra of multiple systems 
can also derive t h e Dop pler shifts, such as the method by 
ISimon k Sturml (119941) that works in wavelength space, 
and that of iHadraval ijl995h in Fourier space, although 
we have not yet seen an application to quadruple-lined 
spectra. 

An implicit requirement of the method developed in 
this paper is that the four templates present a close 
match to the individual spectra of the components. In 
our experience with the use of synthetic templates the 
most important parameters are the effective tempera- 
ture and the rotational broadening. Template mismatch 
is a well known issue in digital cross-correlation that 
can result in bias es in the radial velocities (see, e.g., 
iGriffin et al.1 12000). Similar concerns hold in the four- 
dimensional case (as in the two- and three-dimensional 
versions of TODCOR). Therefore, selection of the op- 
timal templates must be done with care. It is pos- 
sible in many cases to determine the template pa- 
rameters from the observed spectra themselves, if the 
SNR is sufficient. Illustrations of this for the two- 
and th ree-d imensional ve r sions are given by iTorres et al.l 
( 2002) and ITorres et all (|2006ft . The present example 
of HD 110555 does not lend itself to these determina- 
tions because of the faint fourth component, which is 



why we resorted to the use of models. Another possibil- 
ity is to combine our method with one of the disentan- 
gling techniques, which should be able to produce sep- 
arate spectra for each of the components that will have 
higher SNRs than any of the individual spectra. These 
could then be used as templates for the four-dimensional 
cross-correlation scheme in order to derive the velocities. 
While this would seem to be a rather obvious and elegant 
approach, it requires testing to determine the sensitivity 
of the velocities to the details of the disentangling. 
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APPENDIX 



The extension of TO DCO R to four dimensio ns parallels the mathematical development in the Appendices of 
IZucker fc Mazehl (|1994l ) and IZucker et al.l ([1995) for the two- and three-dimensional cases. The reader is referred 
to those papers for some of the details. The starting point is the definition of the one-dimension al cross-correlation 
function between the object spectrum f(n) and the template g(n) (see, e.g.. lTonrv fc Davislll979t ). 

NCTfCTg 

in which N is the number of bins in the observed spectrum, a f and a g are the intensity standard deviations of the 
object and template spectra given by 



N ^— ' x ' ' 3 N 

and s is the shift between the two. As is well known, the use of the Fast Fourier Transform (FFT) allows for a very 
efficient calculation of the numerator of eq.([T]), which is what makes the cross-correlation method practical. This same 
property is key to the success of TODCOR and its extensions. 

We now consider the template g to be a linear combination of four separate templates, one for each star, with 
coefficients a, (3, and 7 being the light ratios of the secondary, tertiary, and quaternary components relative to the 
primary: 

g = gi(n - si) + a g 2 {n - s 2 ) + f3 g 3 (n - s 3 ) + 7 g 4 (n - s 4 ) . 

Substitution in eq.JTJ) yields 

n {A) = E f(n) [gi(n - s^+a g 2 (n - s 2 ) + (3 g 3 {n - s 3 ) + 7 ,g 4 (n - s 4 )] 

Na f a g (s 1 ,s 2 ,s 3 ,s 4 ) 

in which 

S2,s 3 ,s 4 ) = -rp y,[gi(n -sx)+a g 2 (n - s 2 ) + (3 g 3 (n - s 3 ) + 7 g 4 (n - s 4 )] 2 . 
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The numerator of eq.((2]) can be computed efficiently using FFT, since the four terms have the same form as eq.(fT]). 
The standard deviation cr s (si, s 2 , S3, S4) in the denominator is 

o* = a 2 gi + a 2 a 2 2 + /3 2 <7 2 g3 + j 2 a 2 gi + 2aa 12 + 2/3a 13 + 2 1 a 14 + 2a0a 23 + 2a 7 cr 2 4 + 20ja 34 , (3) 

where we have defined, as in IZucker et all (fl995h . 

aij {sj - s. t ) = jj^2gi(n- Si)9j (n - Sj) . 

The first four terms on the right-hand side of eq.([3]) include the standard deviations of the templates. The other six 
terms have again the same form as the numerator in eq.([T]), so they can be computed easily using FFT. If we now 
adopt the following additional definitions 

Ci( Si ) = — ^2f(n) gi (n-Si) (4) 

1 



we may write 



or 



C y (sj - Si) = -—^ V gi (n)gj [n - (sj - si}] , (5) 

N<J g 1 < T g ] 

a l =< + a2(J 92 + P 2a l + 7 2 < + 2(aa gi a g2 C 12 + f3a gi a g3 C 13 + 
7 cr gi cr 94Ci4 + a>f3cr g2 a g3 C 23 + a^a g2 a g4 C 24 + filer g3 a gi C 34 ) , 

o 2 g = * 2 gi [1 + a' 2 + 0' 2 + i 2 + 2(a'C 12 + (3'C 13 + j'C u + a'(3'C 23 + a'j'C 24 + 0'j'Cu)] , 

where a' = cn(a g2 /a gi ), 0' = 0(a g3 /a gi ), and 7' = r y(cr gi /a gi ). After some algebra the final expression for the 
four-dimensional CCF becomes 

n (4) = C 1 +a'C2+0'C 3 +'fC 4 , 

y/l + a' 2 + 0' 2 + 7 ' 2 + 2{a'C 12 + /3>C 13 + 7 'C 14 + a>f3'C 23 + a'YC 24 + 0'YC 34 ) ' 

Thus, for known values of the three light ratios, IZ^ = TZ^(si, s 2 , S3, S4) is reduced to a combination of ten onc- 
dimensional functions that can be easily computed (eq.[4] and eq. [5]): four of them (Cj, i = 1, ...,4) are of the 
observed spectrum against each of the templates, and the other six (Cy , i 7^ j) are the pairwise correlations of the 
four templates against each other. 

In most cases the light ratios of the components are not known a priori, but in fact they can be of considerable 
interest in themselves, as illustrated by the exampl e described in this p aper. It is therefore important to be able to 
estimate them from the observed spectra. Following I Zucker et alj (|1995t ) the light ratios may be computed directly by 
seeking the maximum of the correlation function IZ^ at each combination of shifts {s\, s 2 , s 3 , s^}. This results in a 
system of three equations (dlZ^/da = 0, dTZ^/80 = 0, dlZ^/d^ = 0) with three unknowns (a, 0, 7) which, after 
considerable manipulation, can be solved analytically. We obtain the expressions 

T 9i 



C X A 2 - 


f- c 2 e 4 - 


f- C 3 e 4 - 


hc 4 e 2 


CiAiH 


- C 2 A 2 - 


h C3A3 - 


1-C4A4 


C1A3- 


f- C 2 9i - 


1- c 3 e 5 - 


hc 4 e 3 


CiAiH 


- C 2 A 2 - 


h C3A3 - 


1-C4A4 


C1A4- 


f c 2 e 2 - 


1- c 3 e 3 - 


hc 4 e 6 


CiAiH 


- C 2 A 2 - 


h C3A3 - 


1-C4A4 



/?= p» * . 1 li A ° * * (7) 



where we have defined 



e x 


= C\ 2 C\ 3 


+ C\ 4 C 23 + C 24 C 34 — C12C14C34 


— C\ 3 C\ 4 C 24 


— c 23 


e 2 


= C\ 2 G\4 


+ C\ 3 C 24 + C 23 C 34 — C\ 2 C\ 3 C 34 


— C\ 3 C\ 4 C 23 


— c 24 


e 3 


= C13C14 


+ C 2 2 C 3 4 + C 2 3C 2 4 — C12C13C24 


— C\ 2 C\ 4 G 23 


— C34 


e 4 


= l-C 2 13 


— C 2 4 — C 34 + 


2C13C14C34 






e 5 


= 1-Cl 2 2 


— C 14 — C 24 + 


2Ci 2 Ci4C 2 4 






e 6 


= i-c 2 12 


— Cl3 — ^23 + 


2Ci 2 Ci3C 2 3 







and 



Ai = 1 — C 23 — C24 — C34 + 2C 2 3C 2 4(734 

A 2 = C12C34 + Ci3C 2 3 + Ci4C 2 4 — C13CMC34 — Ci4C 2 3C34 — C\ 2 
A3 = Ci 2 C 2 3 + C13C24 + C14C34 — Ci 2 C 2 4C34 — Ci4C 2 3C 2 4 — C\ 3 
A4 = C\ 2 C 24 + C13C34 + C\ 4 C\ 3 — Ci 2 C 2 3C34 — Ci3C 2 3C 2 4 — C14 
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After rescaling to convert to a', /?', and 7', as indicated above, the values obtained from eqs.0 can be substituted 
in eq. ((6]) to compute the value of the correlation for the optimal values of the three light ratios at each set of shifts 
{si, s 2 , S3, S4}. 
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TABLE 1 

Tests of the algorithm using synthetic spectra. 



Component 


Radial velocity 
(kms -1 ) 


Light ratio 

(a, 0, 7) 


Input values 






Primary 

Secondary 

Tertiary 

Quaternary 


+20.00 
-20.00 
+10.00 
-5.00 


0.50 
0.80 
0.50 


Output from Test 1: 


SNR = 25 




Primary 

Secondary 

Tertiary 

Quaternary 


+20.02 ± 0.47 
-19.67 ± 0.38 
+10.13 ± 0.86 
-4.61 ± 0.75 


0.49 ± 0.06 
0.82 ± 0.14 
0.50 ± 0.06 


Output from Test 2: SNR = 10 




Primary 

Secondary 

Tertiary 

Quaternary 


+20.17 ± 1.44 
-19.37 ± 1.12 
+9.97 ± 2.54 
-4.69 ± 2.24 


0.56 ± 0.27 
1.01 ± 0.70 
0.57 ± 0.32 



TABLE 2 

Results of a model fit to the observational constraints for HD 110555, 

IN ORDER TO INFER THE EFFECTIVE TEMPERATURES OF THE FOUR COMPONENTS. 



Parameter Observed value Model result {0 — C)/cr 



a 


0.92 


± 


0.03 


0.90 a 


+0.67 


P 


0.38 


± 


0.03 


0.41 a 


-1.00 


7 


0.06 


± 


0.03 


0.10 a 


-1.33 


Combined V (mag) 


8.36 


± 


0.01 b 


8.358 


+0.20 


Combined B — V (mag) 


0.658 


± 


0.019 b 


0.667 


-0.47 


Combined V — J (mag) 


1.157 


± 


0.029 c 


1.160 


-0.10 


Combined V — H (mag) 


1.525 


± 


0.037 c 


1.514 


+0.30 


Combined V — K (mag) 


1.552 


± 


0.031 c 


1.552 


0.00 


AV (mag) 


1.33 


± 


0.04 d 


1.383 


-1.32 


7r (mas) 


12.45 


± 


1.82 b 


11.05 


+0.77 



Note. — The m odel fit is based on a 3 Gyr isochrone from the series by 
Girardi ct al. (2000) for solar metallicity. The mass ratios for the two spectro- 
scopic binaries have been adopted from Table [6] as additional constraints. As an 
additional check, the parallax resulting from the models (average of estimates in 
four passbands) is inferred by comparison of the predicted integrated absolute 
magnitudes in VJHK with the apparent magnitudes. 

a Values interpolated between the B and V bands to match the mean wavelength 
of the obse rved light ratios (5188.5 A). b As listed in the Hipparcos Catalogue 
ESA| |T997). C Derive d from the John s on V magnitude and JHK S from 2MASS, 
and converted to the Bcsscll & Brett (1988) photometric system adopted for the 
isochrones using the transformations by Carpenter (2001). d The original Hippar- 
cos measurement of the magnitude difference between the visual components of 
HD 110555 in the H p band (AH V = 1.35 + 0.03) h as been transformed here to 
the V band using the relations by Harmancc (1998). 



TABLE 3 

Properties for the four stars in HD 110555 

INFERRED FROM MODELS. 





Mass 


T e g 


log 3 


My 


Component 


(M ) 


(K) 


(cgs) 


(mag) 


Primary (Aa) 


1.065 


5944 


4.393 


4.540 


Secondary (Ab) 


1.048 


5891 


4.414 


4.655 


Tertiary (Ba) 


0.930 


5457 


4.527 


5.475 


Quaternary (Bb) 


0.758 


4686 


4.634 


6.947 



Note. — Results are based on a comparison of the ob- 
served magnitu des and colors of HD 110555 with a model 
isochrone from Girardi ct al. (2000), for solar composition 
and a representative age of 3 Gyr (see text and Table [2]l . 
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TABLE 4 

Radial velocity measurements for the visual primary (A = Aa + Ab) of HD 110555, in the 

HELIOCENTRIC FRAME. 



HJD RV Aa TAa (0-C) Aa RV Ab a Ab (0-C) Ab 

(2,400,000+) (kms -1 ) (kms -1 ) (kms -1 ) (kms -1 ) (kms -1 ) (kins- 1 ) Phase 



48997.8913 -17.42 1.62 -1.80 +14.99 1.66 +0.86 0.8844 

49025.7939 +0.91 1.22 +0.85 -2.60 1.25 -0.80 0.3711 

49046.7522 -30.05 1.56 -0.28 +28.76 1.60 +0.24 0.7368 

49058.8167 +5.61 1.69 +0.73 -8.37 1.73 -1.66 0.9472 

49062.8294 +31.03 1.46 -0.16 -30.83 1.50 +2.62 0.0172 



Note. — Table [4] is published in its entirety in the electronic edition of the Astrophysical Journal. 
A portion is shown here for guidance regarding its form and content. 



TABLE 5 

Radial velocity measurements for the visual secondary (B = Ba + Bb) of HD 110555, in 

THE HELIOCENTRIC FRAME. 



HJD 


RV B a 


°"Ba 


(O - C) B a 


RV Bb 


o" Bb 


(O - C) Bb 




(2,400,000+) 


(kms 1 ) 


(kms -1 ) 


(kms 1 ) 


(kms 1 ) 


(kms -1 ) 


(kms 1 ) 


Phase 


48997.8913 


-16.59 


3.81 


+0.88 


+14.87 


12.16 


-4.05 


0.3364 


49025.7939 


+13.21 


2.86 


+4.37 


-17.32 


9.14 


-3.96 


0.6999 


49046.7522 


+24.64 


3.67 


-4.86 


-48.85 


11.71 


-10.14 


0.9730 


49058.8167 


-24.71 


3.97 


+0.86 


+22.26 


12.65 


-6.61 


0.1302 


49062.8294 


-28.42 


3.43 


-3.15 


+39.04 


10.96 


+10.54 


0.1825 



Note. — Table IB1 is published in its entirety in the electronic edition of the Astrophysical Journal. 
A portion is shown here for guidance regarding its form and content. 



TABLE 6 

Orbital elements for the two visual components of HD 110555. 



Parameter 



Visual Primary (Aa+Ab) Visual Secondary (Ba+Bb) 



Adjusted quantities 

P (days) 

V (kms" 1 ) 

ifnrim (klllS -1 ). . 



l pnm (kms 

K scc (kms -1 ) . 



"Vim (deg) 

T (HJD-2,400,000) 

Derived quantities 



Mprmi Sin 

Msec sin 3 
: Mi 



3,' 



(M©) 

i (M ) 

q = M BCC /M pliln 

QprimSini (10 6 km) 

dsecsini (10 6 km) 

a sini (Rq ) 

Other quantities pertaining to the fit 

Time span (days) 

Orbital cycles 

Oprim (kms -1 ) 

<t scc (kms -1 ) 



57.32244 


± 


0.00095 


76.7489 


± 


0.0059 


-0.865 


± 


0.069 


-1.13 


± 


0.22 


35.07 


± 


0.15 


30.47 


± 


0.40 


35.66 


± 


0.15 


37.39 


± 


1.17 


0.3030 


± 


0.0028 


0.4970 
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Fig. 1. — Test of the algorithm on an artificial quadruple-lined spectrum with strongly blended lines, (a) One-dimensional CCF of 
a synthetic composite spectrum (GOV + G2V + GOV + G2V) against the template corresponding to the primary star (see text). The 
shifts imposed on the primary, secondary, tertiary, and quaternary are indicated with the vertical dotted lines, and are not measurable 
with standard techniques due to severe blending of the correlation peaks, (b) Cross-section of the four-dimensional CCF taken at its 
maximum, as a function of the primary velocity, with the velocities of the other three components held fixed at the values that maximize 
the correlation, (c)-(e) Same as above, for the other three components. 
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Fig. 2. — Subset of the observations of HD 110555 showing a few of our spectra (left) and the corresponding one-dimensional CCFs 
(right). Dates of observation are as labeled. Occasionally three peaks are seen in the correlation functions, although more commonly only 
one or two are visible. 
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Fig. 3. — (a) One-dimensional cross-correlation of one of our spectra of HD 110555 (shown at the bottom of Fig. [2]l against the 
template corresponding to the primary star. The velocities predicted from the spectroscopic orbits for the primary, secondary, tertiary, and 
quaternary are indicated with the vertical dotted lines, (b) Cross-section of the four-dimensional CCF taken at its maximum, as a function 
of the primary velocity, with the velocities of the other three components held fixed at the values that maximize the correlation, (c)— (e) 
Same as above, for the other three components. 
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Fig. 4. — Schematic location of the components of HD 110555 in the H-R diagram according to our modeling, shown along with an 
isochrone from the model series by Girardi ct al. (2000) for solar composition and a representative age of 3 Gyr. 
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Fig. 5. — Radial velocity measurements and computed curves for each of the visual components of HD 110555. The primary star in 
each spectroscopic binary is represented with filled circles. The centcr-of-mass velocities are indicated with dotted lines. Error bars for the 
visual primary components are typically smaller than the size of the points. 



